######################################################################
##                                                                  ## 
##      PROJECT TITLE:    Informals dataset paper                   ##
##      PROJECT AUTHOR:   Sam S. Rowan; Charles B. Roger            ##
##      EMAIL:            SAM.ROWAN@CONCORDIA.CA                    ##
##                        CROGER@IBEI.ORG                           ##
##      DESCRIPTION:      Code for descriptive statistics and plots ##
##                          for informals dataset paper             ##
######################################################################

#### Init ####
library(haven)
library(tidyverse)
library(directlabels)
library(countrycode)

### Set working directory to parent directory (above 'code')

#####################################################################

# Code to create:
#' Table: Top/bottom IO membership
#' Figure: IO->state membership distribution
#' Figure: State->IO membership distribution 
#' Figure: Map of informals as a percentage of total membership
#' Figure: Modularity over time
#' Table: Summary statistics of IO membership

#####################################################################

## Load in and clean data

# Load in formal IO data 
formal_membership <- read_dta("data/formals_v3.dta") %>%
  filter(year > 1954, year < 2011) %>%
  mutate(iso3c = countrycode::countrycode(.$ccode, 
                                          origin = 'cown',
                                          destination = 'iso3c')) %>% 
  filter(!is.na(iso3c)) %>% 
  pivot_longer(names_to = "io", 
               values_to = "member",
               aaaid:wtouro) %>% 
  select(iso3c, year, io, member) %>% 
  mutate(member = ifelse(.$member<0, NA, 
                         ifelse(.$member>=2, 0, .$member))) %>% 
  group_by(iso3c, year) %>% 
  mutate(formals_sum = sum(member, na.rm = T)) %>% 
  select(iso3c, year, formals_sum) %>% 
  distinct() %>% 
  as.data.frame()

# Load in informal IO data 
informal_membership <- read_dta("data/informals_v3.dta") %>%
  filter(year > 1954, year < 2011) %>%
  mutate(iso3c = countrycode::countrycode(.$ccode, 
                                          origin = 'cown',
                                          destination = 'iso3c')) %>% 
  filter(!is.na(iso3c)) %>% 
  pivot_longer(names_to = "io", 
               values_to = "member",
               acd:zc) %>% 
  select(iso3c, year, io, member) %>% 
  mutate(member = ifelse(.$member<0, NA, 
                         ifelse(.$member>=2, 0, .$member))) %>% 
  group_by(iso3c, year) %>% 
  mutate(informals_sum = sum(member, na.rm = T)) %>% 
  select(iso3c, year, informals_sum) %>% 
  distinct() %>% 
  as.data.frame()

# Combine IO data and make total measure
io_membership <- full_join(formal_membership, informal_membership, 
                           by = c("iso3c", "year")) %>% 
  mutate(total_sum = formals_sum + informals_sum,
         cc_cname = countrycode::countrycode(.$iso3c, 
                                             origin = 'iso3c', 
                                             destination = 'country.name.en'))


#### Tables of the top and bottom joiners, across IO types

## Top joiners

# Total IOs
top_total_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  arrange(-total_sum) %>%
  mutate(rank_total = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_total <= 10) %>%
  arrange(rank_total) %>%
  select(cc_cname, total_sum) %>%
  # select(iso3c, total_sum, rank_total) %>% 
  as.data.frame() 

top_share_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  mutate(informals_share = 100*round((informals_sum/(formals_sum + informals_sum)),3 )) %>% 
  arrange(-informals_share) %>%
  mutate(rank_share = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_share <= 10) %>%
  arrange(rank_share) %>%
  select(cc_cname, informals_share) %>%
  # select(iso3c, total_sum, rank_total) %>% 
  as.data.frame() 
         
         
# Formal IOs
top_formals_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  arrange(-formals_sum) %>%
  mutate(rank_formal = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_formal <= 10) %>%
  arrange(rank_formal) %>%
  select(cc_cname, formals_sum) %>%
  # select(iso3c, formals_sum, rank_formal) %>%
  as.data.frame() 

# Informal IOs
top_informals_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  arrange(-informals_sum) %>%
  mutate(rank_informal = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_informal <= 10) %>%
  arrange(rank_informal) %>%
  select(cc_cname, informals_sum) %>%
  # select(iso3c, informals_sum, rank_informal) %>% 
  as.data.frame() 

# Have a look
top_table <- cbind.data.frame(top_formals_joiners, 
                              top_informals_joiners, 
                              top_share_joiners) %>% 
  magrittr::set_colnames(c("Country", "Formal IOs", 
                           "Country", "Informal IOs", 
                           "Country", "Informals %"))

## Bottom joiners

# Total IOs
bottom_formals_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  arrange(formals_sum) %>%
  mutate(rank_formal = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_formal <= 10) %>%
  arrange(-rank_formal) %>%
  select(cc_cname, formals_sum) %>%
  # select(iso3c, formals_sum, rank_formal) %>% 
  as.data.frame() 

# Formal IOs
bottom_informals_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  arrange(informals_sum) %>%
  mutate(rank_informal = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_informal <= 10) %>%
  arrange(-rank_informal) %>%
  select(cc_cname, informals_sum) %>%
  # select(iso3c, informals_sum, rank_informal) %>% 
  as.data.frame() 

# Informal IOs
bottom_total_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  arrange(total_sum) %>%
  mutate(rank_total = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_total <= 10) %>%
  arrange(-rank_total) %>%
  select(cc_cname, total_sum) %>%
  # select(iso3c, informals_sum, rank_total) %>% 
  as.data.frame() 

bottom_share_joiners <- io_membership %>%
  filter(year == 2010, total_sum!="NA") %>%
  mutate(informals_share = 100*round((informals_sum/(formals_sum + informals_sum)),3 )) %>% 
  arrange(informals_share) %>%
  mutate(rank_share = seq(1, length(unique(.$iso3c)), 1)) %>%
  filter(rank_share <= 10) %>%
  arrange(-rank_share) %>%
  select(cc_cname, informals_share) %>%
  # select(iso3c, total_sum, rank_total) %>% 
  as.data.frame() 


bottom_table <- cbind.data.frame(bottom_formals_joiners, 
                              bottom_informals_joiners, 
                              bottom_share_joiners) %>% 
  magrittr::set_colnames(c("Country", "Formal IOs", 
                           "Country", 
                           "Informal IOs", 
                           "Country", "Informals %"))

## Combine top and bottom into one table

bind_rows(top_table, bottom_table) %>% 
  knitr::kable(format = "html") %>% 
  kableExtra::save_kable(file = "tables/state_io_portfolios.html")
# Then paste this HTML table into Word



#### Quantiles of membership: at the state- and IO-level

## Read in and tidy formal IOs data
formals_gather <- read_dta("data/formals_v3.dta") %>%
  unite("ccode_year", ccode, year, sep = "_") %>%
  gather(key = io_name, value = member, aaaid:wtouro) %>%
  mutate(member = ifelse(member < 0, NA, member)) %>%
  separate(ccode_year, into = c("ccode", "year")) 

# Get quantiles
formals_ts_analysis_iolevel <- formals_gather %>%
  filter(year > 1964, year < 2011, member>0) %>%
  group_by(io_name, year) %>% 
  summarise(sum_io_members = sum(member, na.rm=T)) %>%
  ungroup() %>%
  group_by(year) %>%
  mutate(quantile_10 = quantile(sum_io_members, 0.1, na.rm=T)) %>%
  mutate(quantile_20 = quantile(sum_io_members, 0.2, na.rm=T)) %>%
  mutate(quantile_30 = quantile(sum_io_members, 0.3, na.rm=T)) %>%
  mutate(quantile_40 = quantile(sum_io_members, 0.4, na.rm=T)) %>%
  mutate(quantile_50 = quantile(sum_io_members, 0.5, na.rm=T)) %>%
  mutate(quantile_60 = quantile(sum_io_members, 0.6, na.rm=T)) %>%
  mutate(quantile_70 = quantile(sum_io_members, 0.7, na.rm=T)) %>%
  mutate(quantile_80 = quantile(sum_io_members, 0.8, na.rm=T)) %>%
  mutate(quantile_90 = quantile(sum_io_members, 0.9, na.rm=T)) %>%
  mutate(quantile_max = quantile(sum_io_members, 1, na.rm=T)) %>%
  select(-io_name, -sum_io_members) %>%
  distinct() %>% 
  arrange(year) %>% 
  gather(key = measure, 
         value = members, 
         quantile_10:quantile_max, 
         factor_key = T) %>%
  ungroup() %>%
  mutate(year = as.numeric(year)) %>% 
  as.data.frame() 
  
# Plot formals
formals_ts_analysis_iolevel %>%
  mutate(measure = paste(str_sub(measure, start = -2L, end = -1L),
                         "th", sep = "")) %>% 
  mutate(measure = ifelse(measure == "axth", "Max", measure)) %>% 
  ggplot(aes(year, members, group = measure, color = measure)) + 
  geom_line() +
  scale_color_manual(values = c("grey90", "grey80",
                                "grey70", "grey60",
                                "grey50", "grey40",
                                "grey30", "grey20",
                                "grey10", "black")) +
  coord_cartesian(ylim=c(0,215)) +
  labs(x="Year", y="Number of members", title = "Formal IOs") +
  theme(legend.position = "none", 
                         panel.background = element_rect(fill=NA), 
                         text=element_text(size=18),
                         axis.text = element_text(size = rel(1.5)),
                         panel.grid.major.y = element_line(colour = "grey80"),
                         panel.grid.major.x = element_line(colour = "grey80"),
                         panel.ontop = FALSE,
                         axis.line = element_line(size = .5, colour = "black")) +
  geom_dl(aes(label = measure), method = list(dl.combine("last.points"), cex = 0.8)) +
  geom_point(data = formals_ts_analysis_iolevel %>%
               filter(measure == "quantile_50"),
             aes(year, members),
             shape = 16,
             size = 2,
             color = "black")
ggsave("figures/io_formals_line.pdf", units = 'in', height=6, width=8)

## Read in and tidy informal IOs data
informals_gather <- read_dta("data/informals_v3.dta") %>% 
  select(-state_name, -active) %>%
  unite("ccode_year", ccode, year, sep = "_") %>%
  gather(key = io_name, value = member, acd:zc) %>%
  mutate(member = ifelse(member < 0, NA, member)) %>%
  separate(ccode_year, into = c("ccode", "year")) %>%
  filter(year > 1964, year < 2011) %>%
  arrange(ccode, year, io_name)

## Get quantiles
informals_ts_analysis_iolevel <- informals_gather %>%
  filter(year > 1964, year < 2011, member>0) %>%
  group_by(io_name, year) %>% 
  summarise(sum_io_members = sum(member, na.rm=T)) %>%
  ungroup() %>%
  group_by(year) %>%
  mutate(quantile_10 = quantile(sum_io_members, 0.1, na.rm=T)) %>%
  mutate(quantile_20 = quantile(sum_io_members, 0.2, na.rm=T)) %>%
  mutate(quantile_30 = quantile(sum_io_members, 0.3, na.rm=T)) %>%
  mutate(quantile_40 = quantile(sum_io_members, 0.4, na.rm=T)) %>%
  mutate(quantile_50 = quantile(sum_io_members, 0.5, na.rm=T)) %>%
  mutate(quantile_60 = quantile(sum_io_members, 0.6, na.rm=T)) %>%
  mutate(quantile_70 = quantile(sum_io_members, 0.7, na.rm=T)) %>%
  mutate(quantile_80 = quantile(sum_io_members, 0.8, na.rm=T)) %>%
  mutate(quantile_90 = quantile(sum_io_members, 0.9, na.rm=T)) %>%
  mutate(quantile_max = quantile(sum_io_members, 1, na.rm=T)) %>%
  select(-io_name, -sum_io_members) %>%
  distinct() %>% 
  arrange(year) %>% 
  gather(key = measure, 
         value = members, 
         quantile_10:quantile_max, 
         factor_key = T) %>%
  ungroup() %>%
  mutate(year = as.numeric(year)) %>% 
  as.data.frame() 

informals_ts_analysis_iolevel %>%
  mutate(measure = paste(str_sub(measure, start = -2L, end = -1L),
                         "th", sep = "")) %>% 
  mutate(measure = ifelse(measure == "axth", "Max", measure)) %>% 
  ggplot(aes(year, members, group = measure, color = measure)) + 
  geom_line() + 
  scale_color_manual(values = c("grey90", "grey80",
                                "grey70", "grey60",
                                "grey50", "grey40",
                                "grey30", "grey20",
                                "grey10", "black")) +
  coord_cartesian(ylim=c(0,215)) +
  labs(x="Year", y="Number of members", title="Informal IOs") +
  theme(legend.position = "none", 
        panel.background = element_rect(fill=NA), 
        text=element_text(size=18),
        axis.text = element_text(size = rel(1.5)),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.major.x = element_line(colour = "grey80"),
        panel.ontop = FALSE,
        axis.line = element_line(size = .5, colour = "black")) +
  geom_dl(aes(label = measure), method = list(dl.combine("last.points"), cex = 0.8)) +
  geom_point(data = informals_ts_analysis_iolevel %>%
               filter(measure == "quantile_50"),
             aes(year, members),
             shape = 16,
             size = 2,
             color = "black")
ggsave("figures/io_informals_line.pdf", units = 'in', height=6, width=8)

## State level

formals_ts_analysis_statelevel <- formals_gather %>%
  filter(year > 1964, year < 2011, member>0, ccode!="NA") %>%
  group_by(ccode, year) %>% 
  summarise(sum_io_members = sum(member, na.rm=T)) %>%
  ungroup() %>%
  group_by(year) %>%
  mutate(quantile_10 = quantile(sum_io_members, 0.1, na.rm=T)) %>%
  mutate(quantile_20 = quantile(sum_io_members, 0.2, na.rm=T)) %>%
  mutate(quantile_30 = quantile(sum_io_members, 0.3, na.rm=T)) %>%
  mutate(quantile_40 = quantile(sum_io_members, 0.4, na.rm=T)) %>%
  mutate(quantile_50 = quantile(sum_io_members, 0.5, na.rm=T)) %>%
  mutate(quantile_60 = quantile(sum_io_members, 0.6, na.rm=T)) %>%
  mutate(quantile_70 = quantile(sum_io_members, 0.7, na.rm=T)) %>%
  mutate(quantile_80 = quantile(sum_io_members, 0.8, na.rm=T)) %>%
  mutate(quantile_90 = quantile(sum_io_members, 0.9, na.rm=T)) %>%
  mutate(quantile_max = quantile(sum_io_members, 1, na.rm=T)) %>%
  select(-ccode, -sum_io_members) %>%
  distinct() %>% arrange(year) %>% 
  gather(key = measure, 
         value = members, 
         quantile_10:quantile_max, 
         factor_key = T) %>%
  as.data.frame() %>% 
  ungroup() %>%
  mutate(year = as.numeric(year)) 

formals_ts_analysis_statelevel %>% 
  mutate(measure = paste(str_sub(measure, start = -2L, end = -1L),
                         "th", sep = "")) %>% 
  mutate(measure = ifelse(measure == "axth", "Max", measure)) %>% 
  ggplot(aes(year, members, colour = measure, fill = measure)) +
  geom_line() +
  scale_color_manual(values = c("grey90", "grey80",
                                "grey70", "grey60",
                                "grey50", "grey40",
                                "grey30", "grey20",
                                "grey10", "black")) +
  coord_cartesian(ylim=c(0,125)) +
  labs(x="Year", y="Number of memberships", title="Formal IOs") +
  theme(legend.position = "none", 
        panel.background = element_rect(fill=NA), 
        text=element_text(size=18),
        axis.text = element_text(size = rel(1.5)),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.major.x = element_line(colour = "grey80"),
        panel.ontop = FALSE,
        axis.line = element_line(size = .5, colour = "black")) +
  geom_dl(aes(label = measure), method = list(dl.combine("last.points"), cex = 0.8)) +
  geom_point(data = formals_ts_analysis_statelevel %>%
               filter(measure == "quantile_50"),
             aes(year, members),
             shape = 16,
             size = 2,
             color = "black")
ggsave("figures/state_formals_line.pdf", units = 'in', height=6, width=8)

informals_ts_analysis_statelevel <- informals_gather %>%
  filter(year > 1964, member>0) %>%
  group_by(ccode, year) %>% 
  summarise(sum_io_members = sum(member, na.rm=T)) %>%
  ungroup() %>%
  group_by(year) %>%
  mutate(quantile_10 = quantile(sum_io_members, 0.1, na.rm=T)) %>%
  mutate(quantile_20 = quantile(sum_io_members, 0.2, na.rm=T)) %>%
  mutate(quantile_30 = quantile(sum_io_members, 0.3, na.rm=T)) %>%
  mutate(quantile_40 = quantile(sum_io_members, 0.4, na.rm=T)) %>%
  mutate(quantile_50 = quantile(sum_io_members, 0.5, na.rm=T)) %>%
  mutate(quantile_60 = quantile(sum_io_members, 0.6, na.rm=T)) %>%
  mutate(quantile_70 = quantile(sum_io_members, 0.7, na.rm=T)) %>%
  mutate(quantile_80 = quantile(sum_io_members, 0.8, na.rm=T)) %>%
  mutate(quantile_90 = quantile(sum_io_members, 0.9, na.rm=T)) %>%
  mutate(quantile_max = quantile(sum_io_members, 1, na.rm=T)) %>%
  select(-ccode, -sum_io_members) %>%
  distinct() %>% arrange(year) %>% 
  gather(key = measure, 
         value = members, 
         quantile_10:quantile_max, 
         factor_key = T) %>%
  as.data.frame() %>% 
  ungroup() %>%
  mutate(year = as.numeric(year)) 

informals_ts_analysis_statelevel %>%
  mutate(measure = paste(str_sub(measure, start = -2L, end = -1L),
                         "th", sep = "")) %>% 
  mutate(measure = ifelse(measure == "axth", "Max", measure)) %>% 
  ggplot(aes(year, members, colour = measure, fill = measure)) +
  geom_line() +
  scale_color_manual(values = c("grey90", "grey80",
                                "grey70", "grey60",
                                "grey50", "grey40",
                                "grey30", "grey20",
                                "grey10", "black")) +
  coord_cartesian(ylim=c(0,125)) +
  labs(x="Year", y="Number of memberships", title = "Informal IOs") +
  theme(legend.position = "none", 
        panel.background = element_rect(fill=NA), 
        text=element_text(size=18),
        axis.text = element_text(size = rel(1.5)),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.major.x = element_line(colour = "grey80"),
        panel.ontop = FALSE,
        axis.line = element_line(size = .5, colour = "black")) +
  geom_dl(aes(label = measure), method = list(dl.combine("last.points"), cex = 0.8)) +
  geom_point(data = informals_ts_analysis_statelevel %>%
               filter(measure == "quantile_50"),
             aes(year, members),
             shape = 16,
             size = 2,
             color = "black")
ggsave("figures/state_informals_line.pdf", units = 'in', height=6, width=8) 

#### Modularity ####

modularity <- read_csv("data/ml_summaries_cow_v3.csv")

modularity %>%
  pivot_longer(names_to = "IOs",
               values_to = "modularity",
               ml_formals:ml_total) %>%
  filter(year > 1964) %>%
  mutate(IOs = ifelse(.$IOs == "ml_formals", "Formals",
                      ifelse(.$IOs == "ml_informals", "Informals",
                             "Total"))) %>%
  ggplot(., aes(year, modularity, colour = IOs)) +
  coord_cartesian(ylim=c(0,.25)) +
  labs(x="Year", y="Modularity", title = "") +
  theme(legend.position = "bottom", 
        panel.background = element_rect(fill=NA), 
        text=element_text(size=18),
        axis.text = element_text(size = rel(1.5)),
        panel.grid.major.y = element_line(colour = "grey80"),
        panel.grid.major.x = element_line(colour = "grey80"),
        panel.ontop = FALSE,
        axis.line = element_line(size = .5, colour = "black")) +
  # geom_smooth(method = "loess") +
  geom_line(size = 1) +
  scale_color_manual(values = c("grey60", "grey30", "black")) +
  scale_shape_manual(values=c(3, 16, 17)) +
  geom_point(aes(shape = IOs), size = 3) 
ggsave("figures/io_modularities.pdf", units = 'in', height=6, width=8)

# ----- 3. Map

WorldData <- map_data('world')
WorldData <- WorldData %>% 
  filter(region != "Antarctica") %>% 
  filter(region != "Greenland") 
WorldData <- fortify(WorldData)
WorldData$iso3c <- countrycode(WorldData$region, origin = 'country.name', destination = 'iso3c')


map_df <- io_membership %>% 
  filter(year == 2010) %>% 
  mutate(informal_share = informals_sum/total_sum) %>% 
  select(iso3c, informal_share) %>% 
  full_join(., WorldData, by = "iso3c") %>% 
  arrange(order)

ggplot() +
  geom_map(data=WorldData, map=WorldData,
           aes(x=long, y=lat, group=group, map_id=region),
           fill=NA, colour="#7f7f7f", size=0) + 
  geom_map(data=map_df, map=WorldData,
           aes(fill=informal_share, map_id=region),
           colour="#7f7f7f", size=0) +
  coord_map("rectangular", lat0=0, xlim=c(-180,180), ylim=c(-60, 90)) + 
  scale_fill_continuous(low = "grey90", high = "black", guide="colorbar") +
  scale_y_continuous(breaks=c()) +
  scale_x_continuous(breaks=c()) +
  labs(fill="Informals as % of total IO membership", title="", x="", y="") +
  theme_bw() +
  guides(fill = guide_colourbar(barwidth = 10)) +
  theme(panel.border = element_blank(),
        text = element_text(size=18),
        legend.position = "bottom") #
ggsave("figures/map_share_2010.pdf", units = "in", height=6, width=8)




#### Summary statistics ####

years <- seq(1925, 2010, 5)

formal_summary <- read_dta("data/formals_v3.dta") %>%
  # filter(active == 1) %>%
  mutate(iso3c = countrycode::countrycode(.$ccode, 
                                          origin = 'cown',
                                          destination = 'iso3c')) %>% 
  filter(!is.na(iso3c)) %>% 
  pivot_longer(names_to = "io", 
               values_to = "member",
               aaaid:wtouro) %>% 
  select(iso3c, year, io, member) %>% 
  mutate(member = ifelse(member<0, NA, 
                         ifelse(member>=2, 0, member))) %>% 
  group_by(iso3c, year) %>% 
  mutate(formals_sum = sum(member, na.rm = T)) %>% 
  select(iso3c, year, formals_sum) %>% 
  distinct() %>% 
  filter(year %in% years) %>% 
  group_by(year) %>% 
  mutate(n = length(unique(iso3c)),
         mean = mean(formals_sum, na.rm=T),
         sd = sd(formals_sum, na.rm = T),
         min = min(formals_sum, na.rm = T),
         quantile25 = quantile(formals_sum, 0.25, na.rm=T),
         median = quantile(formals_sum, 0.5, na.rm=T),
         quantile75 = quantile(formals_sum, 0.75, na.rm=T),
         max = max(formals_sum, na.rm = T)) %>% 
  mutate_at(vars(mean:max), ~round(., 1)) %>% 
  select(-iso3c, -formals_sum) %>% 
  mutate(io = "Formal IOs") %>% 
  distinct() 

informal_summary <- read_dta("data/informals_v3.dta") %>%
  filter(active == 1) %>% 
  mutate(iso3c = countrycode::countrycode(.$ccode, 
                                          origin = 'cown',
                                          destination = 'iso3c')) %>% 
  filter(!is.na(iso3c)) %>% 
  pivot_longer(names_to = "io", 
               values_to = "member",
               acd:zc) %>% 
  select(iso3c, year, io, member) %>% 
  mutate(member = ifelse(.$member<0, NA, 
                         ifelse(.$member>=2, 0, .$member))) %>% 
  group_by(iso3c, year) %>% 
  mutate(informals_sum = sum(member, na.rm = T)) %>% 
  select(iso3c, year, informals_sum) %>% 
  distinct() %>% 
  filter(year %in% years) %>% 
  group_by(year) %>% 
  mutate(n = length(unique(iso3c)),
         mean = mean(informals_sum, na.rm=T),
         sd = sd(informals_sum, na.rm = T),
         min = min(informals_sum, na.rm = T),
         quantile25 = quantile(informals_sum, 0.25, na.rm=T),
         median = quantile(informals_sum, 0.5, na.rm=T),
         quantile75 = quantile(informals_sum, 0.75, na.rm=T),
         max = max(informals_sum, na.rm = T)) %>% 
  mutate_at(vars(mean:max), ~round(., 1)) %>% 
  select(-iso3c, -informals_sum) %>% 
  mutate(io = "Informal IOs") %>% 
  distinct() 

rbind(informal_summary,
      formal_summary) %>% 
  select(Year = year, `Num. states` = n, 
         Mean = mean, `Std. Dev.` = sd, 
         `Min.` = min, `25th pctl` = quantile25,
         Median = median, `75th pctl` = quantile75,
         `Max.` = max, IO = io) %>% 
  # knitr::kable(.,
  #              format = "latex",
  #              linesep = "",
  #              booktabs = T,
  #              caption = "Summary statistics")
  knitr::kable(format = "html") %>% 
  kableExtra::save_kable(file = "tables/summary_statistics.html")
# Then paste this HTML table into Word
